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ABSTRACT 

We study the evolution of the phase-space of collisionless N-body systems under re- 
peated stirrings or pertu rbations. We find co nvergence towards a limited solution 
group, in accordance withlHansen et al. (20101), that is independent of the initial sys- 



tem and environmental conditions, paying particular attention to the assumed gravi- 
tational paradigm (Newtonian and MOND). We examine the effects of changes to the 
perturbation scheme and in doing so identify a large group of perturbations featuring 
radial orbit instability (ROI) which always lead to convergence. The attractor is thus 
found to be a robust and reproducible effect under a variety of circumstances. 

Key words: galaxies: haloes, galaxies: kinematics and dynamics, methods: N-body 
simulations, methods: numerical 



1 INTRODUCTION 

One of the ongoing problems in the field of galactic dynam- 
ics is the non-keplerian nature of rotation curves in spiral 



galaxies (Salucci et al. 2007) and, consequently, the inferred 
presence of massive but undetectable structures that enve- 
lope the luminous component of galaxies. According to the 
prevailing cosmological models, this structure is comprised 
of weakly-interacting particles which are known as Dark 
Matter (DM). Of particular interest in the study of DM 
is the characteristic density profiles of the halos that sur- 
round observed baryonic structures as, without insight into 
the forms taken by these halos, it is hard to make mean- 
ingful predictions in the theory. Cosmological N-body sim- 
ulations suggested early on that all stable halos should look 
like isothermal spheres that could be fit by some universal 
profile ( Dubinski &; Carlberg||199lj and references therein). 
That profile, named the 'NFW profile' after Nav arro et al.| 
was found to be a two-power model with parameters 
a = 1, /3 = 3 as follows: 
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Due to both the simplicity of the model and the ap- 
parent universality of the result, NFW profiles have become 
accepted as the 'go-to' model for simulating halo charac- 
teristics with many results being based on them. Because 
of this, it is vital that the ubiquity of the NFW profile be 
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properly understood and the fact that it is not is a cause 
for concern. This is especially true since the universality of 
the profile has been called into question in the past by X- 
ray observations ( Makino et al.||1998 ) and the Tully-Fisher 
relation ( |McGaugn fe De Blok|1998| ). 

One potential problem is that in the original study, 



Navarro et al. (1996|, only fully-equilibriated halos were se- 
lected from the initial, low-resolution simulation for further, 
high-resolution simulation. Subsequent investigations into 
the larger population have shown that while the profile is 
still a decent fit, even for non-equilibrium structures (Jing 



2000), there is evidence that a less global profiling system 



may have to be used instead ( Host fc Hansen|201l| . 

We focus on recent work by |Hansen et al.| ( |2010| ) (here- 
after HJS) where it was suggested that all collisionless sys- 
tems will tend to move towards characteristics drawn from 
narrow range if they are gently perturbed from their cur- 
rent equilibrium and are allowed to find a new one. HJS 
used a simple algorithm to disturb a set of relaxed systems 
multiple times and observed, in each, a tendency for each 
successive equilibrium to converge to a particular region in 
the allowed parameter space. We attempt to test the repro- 
ducibility of this phenomena, examine how easily the be- 
haviour can be disturbed and attempt to quantify and ex- 
plain this behaviour. We do not focus on other cosmological 
simulation work in this paper as our aim here is to repro- 
duce and understand the mechanism behind the attractor. 
It is not the aim of this paper to discuss the appearance or 
otherwise of the attractor in either observation or simula- 
tion as it is premature to discuss the physical relevance of 
an effect that is not understood. Future work will use this 
paper as a foundation to discuss the physical relevance of 
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the attractor after the mechanism and phenomenology are 
fully understood. 



ATTRACTORS IN THE JEANS' 
PARAMETER SPACE 



The Jeans equations (Jeans 1915) describe the relationship 
between a density field, u(x,t), a potential, <E>(x,£) and the 
arrangement of velocity vectors in the system. In most cases 
there is not enough information to find a single, unique so- 
lution for an unknown component and one must accept a 
range of permitted solutions instead. For a system with no 
net velocity, such as the ones used in this work, a useful form 
of the equation is the following: 
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2 is the radial velocity dispersion, a 2 is the tan- 
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gential velocity dispersion 
other terms are as follows: 
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where M tot is the total mass of the system. There are 
limitless configurations for which this equation will hold. It 
is this that makes any convergent result so surprising; that 
in an effectively infinite parameter space only one subset of 
results should be favoured. The result that HJS found was 
a strong link in the parameter space between the quantities 
f3 and 7 + k to which they fit an empirical relationship thus: 
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(4) 



which removes a degree of freedom from the system. 
HJS find that equation [4] is an 'attractor' i.e. a solution that 
systems will converge towards if they are free to move in 
it's parameter space. In order to observe evolution towards 
the attractor we need to perturb an initially equilibriated 
system in some physical meaningful way and see if it settles 
into equilibria successively closer to one particular solution. 
The attractor can be represented as a linear relation when 
projected into the right parameter space such as in figure [I] 
The data in this plot is taken from the final states of the 
simulations in HJS and it is that data which will provide us 
with our attractor curve throughout this work. 

However, since this space is somewhat unintuitive, we 
shall usually plot the attractor in the more accessible pa- 
rameter of space /3 versus 7 + k. In this space the attractor 
has a more complex 'S-bend' shape which shall be seen in 
subsequent plots. 



3 NUMERICAL SIMULATIONS 

3.1 Initial conditions and numerical code 

The process of perturbation and relaxation is carried out 
using N-body methods designed to replicate those used by 
HJS. Initial conditions are generated using the method out- 




Figure 1. 1-D form of the attractor in the parameter space 
plotting HJS's data for /3 against the / (7, k) given in equation [4] 

the distribution function of the system, /(e, L), into two sep- 
arate functions: an energy distribution function, g(e), which 
controls how energy levels are populated and a circularity 
function, h(x) = h( Lq+l L . — ^y), which controls how circu- 
lar the orbits of energy e are. From these two functions one 
can directly produce functions for density and velocity dis- 
persion necessary for building a Jeans'-stable system, for 
example: 
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lined in Gerhard (1991). This approach is based on splitting 



(x§(e,r, L ) - x 2 ) 
(5) 

with comparable formulae for other velocity dispersions. 
Accordingly, a desired profile can be analytically generated 
by the appropriate choice h(x) (with g(e) being subsequently 
derived using Lucy's method) and that profile can be ran- 
domly populated by sampling the /(e, L). 

Our initial systems are Plummer spheres consisting of 
750,000 particles with a total mass of 5 x 10 8 M . We 
chose to only investigate one density profile in detail as it 
allows direct comparison of the effects of different pertur- 
bation schemes. Plummer spheres were chosen as they are 
formally unrelated to the NFW profile and are easy to create 
with varying anisotropics. We create a variety of Plummer 
spheres to test various aspects of the attractor. We also use 
one set of initial conditions designed to be, at any radius, 
more strongly anisotropic than the attractor. This system 
has a power law density gradient of r~ 2 throughout and a 
strongly radially anisotropic profile. 

We make use of the N-body code NMODY, a particle 
mesh code capable of implementing MOND, for our numer- 
ical simulation. NMODY uses an iterative scheme alongside 
a standard particle mesh to compute Newtonian forces and 
implements an extension to these for working in MOND as 
follows. Firstly, the code examines the density distribution 
and approximates it as an exact spherical solution, or takes 
the previous state of the system, as the starting point. The 
difference between the assumed density distribution and the 
actual distribution is assessed by assessing the convergence 
between the acceleration field produced by the density and 
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the density implied by the potential. When the two are, 
to within some user-defined accuracy, producing the same 
accelerations then that potential is used. The code uses a 
standard leapfrog scheme to apply timesteps. For a full de- 



scription of the code see jCiotti et al.] ([20Q6) 

The ability to have simulations run in MOND is used 
to examine whether using a different form of Poisson's equa- 
tion has any impact on the attractor. If the attractor is a 
gravitational effect then it would be expected that it would 
be fundamentally linked to the description of gravity being 
used. Thus, we can test whether this is true by altering the 
fundamental equations governing gravity in our simulation 
and attempting to identify any impact that this has on the 
attractor. Since MOND performs just this kind of funda- 
mental alteration to gravity a simulation run in the MOND 
paradigm should reveal whether the attractor arises from 
gravity or not depending on whether or not the convergence 
is different. It is because of this that MOND is useful in 
narrowing down the mechanism that drives convergence. 

Since MOND is used primarily as a tool for altering 
Poisson's equation we will not deal in depth with the theory 
itself or the mathematics behind it. Suffice to say here that 
MOND strengthens gravity compared to Newtonian grav- 
ity but only in regimes of very low acceleration. Thus we 
have our three main paradigms of Newtonian gravity, 'per- 
turbative' or 'weak' MOND where the accelerations are high 
and the MONDian modifications are only a low order per- 
turbation and 'deep' or 'strong' MOND where the bulk of 
the system is in the low acceleration regime and the MOND 
effect is significant. In order to maintain the same density 
profile and total mass between all three scenarios while still 
operating in the relevant acceleration regimes, the scale radii 
of the models are bigger for models with a stronger MOND 
influence: Newtonian simulations use 0.05 kpc, weak MOND 
also uses 0.05 kpc and deep MOND 1.0 kpc. For more details 
please see Appendix [A Bekenstein & Milgrom (1984), Lon- 
drillo & Nipoti (2008) or the introduction of Ciotti et al. 



( 2006 and references therein) . 



3.2 Basic Perturbation 

After choosing which model to use for our initial conditions 
(ICs) we define a simple algorithm that we hope will give 
us our evolution towards an attractor. Our principle algo- 
rithm is taken from that used by HJS but with some minor 
differences such as having fixed-mass bins rather than fixed- 
radius bins and applying different factors to each velocity 
component rather than the same one three times. The bin- 
ning in the simulation allows us to define conservation laws 
for localised groups of particles in the simulations rather 
than only conserving over the simulation as a whole, soften- 
ing the impact of the perturbation. By contrast, the binning 
in the analysis is construct quantities such as velocity dis- 
persion and anisotropy which break down for single parti- 
cles. For the former we choose to define bins as radial shells 
containing a fixed number of particles as it enables good 
statistics for the conservation laws in the outer edges of the 
system where number densities are lower. Other schemes are 
employed to test various components of the attractors be- 
haviour and they will explained in their respective sections. 



• Set up a series of radial bins. We chose to create bins 
defined to contain 5,000 particles. 

• For each particle in each bin we examine each of the 
three orthogonal velocity vectors and multiply each by a 
random number drawn from a uniform distribution centered 
around unity e.g. 1 — 0.25 < f < 1 + 0.25. This is referred to 
as the 'shock' or 'perturbation' and / can be called the 'kick 
scale factor'. As a shorthand / is expressed as ±n, so the 
previous example would be communicated as / = ±0.25. 

• Make a choice about what quantities to conserve in the 
system. Here we shall only deal with energy conservation 
as we currently have no algorithm that can conserve both 
energy and angular momentum without either failing to con- 
verge or introducing biases. For an example of an equivalent 
procedure for angular momentum please refer to Appendix 

m 

• To conserve energy, the energy in the bin is reassessed 
in order to rescale all the velocity components equally to 
enforce conservation. Note that, since we are not moving 
particles around, only the kinetic components before and 
after, Tinit and Tf ina i, need consideration: 



- final 
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• If a particle has been perturbed such that it is no longer 
bound in the system then it is re-randomised and the con- 
servation algorithm is reapplied. 

• Derive a dynamical timescale for the system 
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where we are using the 95* mass percentile as a represen- 
tative distance for the system. For our initial systems this is 
equivalent to approximately 3 scale radii. 

• The system is then left to evolve in an N-body simula- 
tor for 1 dynamical timescales. This 'flow' period allows the 
system to relax and find a new equilibrium. If we were to 
apply another shock too soon then the impact of the second 
shock would be indistinguishable from that of the first. 

• The entire 'shock- flow' cycle is repeated 30 times. 



3.3 Analysis 

Where possible, information about the system is taken di- 
rectly from the 'per-particle' position and velocity file or the 
diagnostic output from the N-body simulation. Local den- 
sity and velocity dispersion are interpolated onto a spherical 
polar grid. Particles have a linear smoothing kernel applied 
to them and the mass contribution at a point of interest are 
summed to find the local density. In the case of the velocity 
dispersion we apply the same method but weight the parti- 
cles by their velocities while applying the smoothing kernel. 
This smooths the small scale noise that otherwise tends to 
dominate the dispersion. We choose points arranged on a 
spherical grid and average over the angular space to find ra- 
dial profiles. Note that the use of a spherical polar grid here 
means that our analysis uses bins of radius in contrast to 
our simulations which use bins of equal mass. 
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Figure 2. Model particle densities projected along orthogonal 
lines of sight into the system. 

4 RECREATING THE HJS RESULTS 
4.1 Outline of simulations 

Various simulations are run over the course testing many 
aspects of the attractor. Table |4.1| summarises, for clarity, 
what has been run and the order in which it will be dis- 
cussed. 

Terms such as 'bimodal kick' are explained in detail in 
their relevant sections. 
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Figure 3. Evolution of the anisotropy profile of the Newtonian, 
isotropic initial conditions at successive (0-black, 5-blue and 30- 
green) 'shock-flow' cycles with the attractor in red. 



majority of systems retain spherical symmetry there is no 
issue in defining v r unambiguously leaving only the question 
of how to treat v t . Since we are explicitly evaluating vg and 
to calculate the anisotropy we are forced to define a of 
anyway. To do so for a non spherical system is perfectly valid 
and will convey all the information it needs to regardless of 
the system's sphericity or preferred axis of rotation. 



4.2 The assumption of spherical symmetry 

Before any in depth analysis can be made it is important 
to learn whether it is reasonable to compress the available 
phase-space, f(r, 0, 0, v r ,ve, v^), according to spherical sym- 
metry, f(r, v r ,Vt) : as we would like to assume from the sym- 
metries of both our IC's and our perturbation algorithm. 
This is important as failure to account for lack of sphericity 
could cause the system to appear more isotropised than be- 
fore (see Appendix [c]). First we must prove that our pertur- 
bation maintains the spatial symmetry of the model which 
we find in figure [2] 

We find that, in general, our algorithms take systems 
with axis ratios (computed by comparing the radii at which, 
as a function of angle, the density profile reaches a value such 
that half the mass of the system is enclosed) of, for example, 
(1.0000,1.0005,1.0005) and, after 25 perturbations, return 
systems with axis ratios of (1.0000,0.9945,1.0167). In rare 
instances where very dispersed systems undergo rapid col- 
lapse statistical variance in the particles positions can lead 
to the development of triaxiality (see the deep MOND simu- 
lations where in only 5 perturbations the axis ratios change 
from (1.0000,0.9991,0.9992) to (1.0000,0.8681,1.3851)). In 
these instances the systems still possess clear radial profiles 
which shall be discussed but detailed discussion the triaxial- 
ity will be avoided since at present it is not of much relevance 
here. 

Now we need to check that the systems do not develop 
a preferred axis of rotation to be sure it is reasonable to 
speak only of radial and tangential velocity. Since the vast 



4.3 Effect of initial anisotropy profiles 

Equation [4] cannot be rearranged as a linear relationship be- 
tween /3 and 7 + k, yet this is the parameter space in which 
the data is most easily visually interpreted. Consequently, 
it is awkward to use equation [4] to plot the attractor in that 
parameter space despite that being easiest for the reader. 
Accordingly, rather than derive the attractor directly from 
equation [4] we plot data from an exemplar converged solu- 
tion from HJS as that is visually cleaner. This is why the 
attractor solution appears as a series of data points rather 
than a smooth, continuous function or allowed region. We 
first look in detail at the simplest of our initial conditions, 
the isotropic Plummer sphere in Newtonian gravity, expect- 
ing our data to converge around the attractor. 

In figure [3] the evolution of the system can be clearly 
seen. The initial condition is isotropic and is therefore a 
vertical line on f3 = with subsequent shocks apparently 
converging on a final state. The movement towards the end 
state is convergent as successive shocks disrupt the previous 
equilibrium less and less until, at a point around 20-30 shock 
cycles, the system reaches a stable state. 

In interpreting plots such as figure [3] a useful rule of 
thumb is that higher values of 7 + k correspond to the inner 
regions of the system and more negative ones to the outer 
edges. This is not an exact relationship due to statistical 
noise and the opposing gradients of p(r) (increasingly nega- 
tive) and of (r) (increasingly positive) , but can still be useful 
as an approximation due to the domination of the 7 term. 
This tells us that the changing anisotropy becomes apparent 
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Table 1. Outline of all simulations performed. 
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Figure 4. Initial (triangles) and final (circles) states for radial 
(blue) and tangential (green) initial conditions. 



Figure 5. Initial and final states for states that lie above the 
attractor 



first in the outer regions of the system before progressively 
spreading inwards. 

Next we consider whether, given that the attractor 
drives an anisotropy gradient, will an initially anisotropic 
model affect the evolution towards the attractor? 

The anisotropy in the initial conditions is clearly visi- 
ble in figure |as a gradient in the initial data points but 
again the convergence is very apparent. This implies that 
the emergence of the attractor is robust to the initial char- 
acter of the halo which is a feature that any explanation of 
the universal profile needs to have. Tangential models take 
longer to converge as the change in anisotropy needed is 
greater. 

Finally we try conditions that are significantly more 
radially anisotropic than the attractor at any given radius. 
These solutions also converge to the attractor just as the 
others do, as shown in figure [5] 



4.4 Effect of gravity theory 

Simulations were run where the initial conditions and sub- 
sequent evolution of the system were carried out under a 
different assumed theory of gravity during the flow phase. If 
the attractor is a phenomenon strongly tied, as one might 
reasonably expect, to gravitational laws then changing the 
form of those laws may have some impact on the evolution 
of the system. 

We test two instances of MONDian dynamics using the 
same perturbations as before; a perturbative case where sys- 
tems of comparable scale to those seen so far are simply 
translated to their stable MOND equivalent with MONDian 
effects only affecting the outer regions and a more extended 
system of the same mass where the entire system is in the 
MONDian regime (see figure |6]). 

The change to a MONDian context seems to have no 
noticeable impact on the evolution of the system as both 
the isotropic and radially anisotropic models, seen in figures 
[7] and [8] respectively, converge. The convergent solution is 
the same as in previous Newtonian simulations and emerges 
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Figure 6. Initial curves of acceleration (solid lines) with radius 
for Newtonian (black), 'perturbative' MOND (green) and deep 
MOND (blue) models compared to their theoretical values (dot- 
ted lines) as set out by appendix [A] Below the red line represents 
the region in which the MOND effects become roughly of order 
unity (a < 3600 (km/s) 2 /kpc) corresponding to the deep MOND 
region. 



Figure 8. The radial perturbative MONDian model. 




Figure 7. The isotropic perturbative MONDian model. 



over the same timescale which is not unexpected since the 
MOND effect in these systems is only slight. 

When we look at the behaviour of the larger, deep 
MOND system in figure we still see evidence of movement 
towards the attract or as shown in figure [9] This plot is not as 
clear as previous ones due to both increased statistical noise 
due to larger characteristic radii and the collapse of the sys- 
tem leading to triaxiality as mentioned in section [472] The 
collapse is interesting in its own right and is present to some 
extent in all systems perturbed in this manner. This will be 
discussed in more detail in section [6] 

We would not expect that the deep MOND regime 
would cause any significantly different behaviour from the 
weak MOND regime due, once more, to the tendency of 
our perturbation method to cause a degree of contraction in 




Figure 9. The deep MONDian model. 
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Figure 10. The collapse of the deep MOND system after only 5 
perturbations (red) compared to the initial mass profile (black). 
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the system to which it is applied. In this case the collapse 
concentrates the mass and most of the mass in the system 
quickly leaves the deep MOND regime as, as figure [To] shows, 
the system becomes more concentrated. 

We see that there is remarkable agreement amongst the 
simulations as they have all converged to a single solution. 
We see clear evidence for convergence to the same solution 
as that found by HJS. This demonstrates that the details 
of the gravitational interactions between particles are not 
especially important to the attractor. This makes sense when 
viewed the same way as the HJS results whereby the flow 
allows for the phase-mixing of the system. One would not 
expect small changes to gravitational coupling to have any 
significant impact on the system's ability to phase-mix and, 
indeed, this is what we find. Now we attempt to disturb the 
convergence by altering the algorithm. 




Shock number 



5 ALTERING THE ALGORITHM 

5.1 Implicit coordinate system 

The perturbation algorithm that we use has dealt with ve- 
locity vectors in a purely cartesian way however we want 
reassurance that the attractor is independent of the coor- 
dinate system. The Newtonian, isotropic system was re-run 
but this time the random scaling of velocity vectors V[ x , v ,z\ 
was carried out on t>[ r ,6>,</>] instead. When the velocities are 
converted back to into Cartesian format the scaling will no 
longer be uniform or in the same limits as before thanks to 
the interconnected and non-linear nature of the transforms 
used e.g.: 

r = xx yy jl zz _^ j. _ r ^ CO s(0) cos(cj))d — sin(0) sin(0)0 
yx 2 + y 2 + z 2 ^ - 

(8) 

We find that the move from Cartesian to spherical co- 
ordinates has no effect on the convergent solution. This is 
in accordance with the findings of HJS where they tested 
different limits on the magnitude of the random factor used 
to scale the velocities. 

5.2 Random scale factors and flow time 

In order to talk quantitatively about the speed of conver- 
gence and to demonstrate that the systems have actually 
reached a convergent solution we need some way of defin- 
ing convergence. Since the evolution of the system is most 
clearly seen in the changing anisotropy profile we choose to 
quantify the convergence by comparing the anisotropy pro- 
files before and after a kick-flow cycle and comparing that 
with the statistical drift that one would expect to see in the 
anisotropy profile of a system that is in equilibrium. This 
is visualised in terms of a plot such as figure [TT] where, at 
every radius, the change in anisotropy between successive 
steps is plotted and then smoothed into contours to show 
the convergence. Figure [TT] is a control plot which shows 
the results of such a plot on the isotropic Newtonian IC's 
when no shocks are applied. Note that since the changes in 
anisotropy are expressed clearest as logarithms the change 
being evaluated is log \/3f — f3i\. 

Figure [TT] shows us that the equilibriated central regions 



Figure 11. The changes in anisotropy in a stable, unperturbed 
system. Note that since the changes in anisotropy are expressed 
clearest as logarithms the change being evaluated is log \(3f — 
The 'shocks' listed on the x-axis are only marked for easy time 
comparison with plots where shocks were actually applied. Here 
the shocks simply mark intervals of 1 T^ yn . Note the difference 
in scale from figure 

are very stable with fluctuations on the scale of ±10 -4 in- 
creasing to dzlO -2 towards the edge of the system. Thus we 
can define convergence in terms of a change in beta that is 
less than or equal to the change at the same radius in the 
unperturbed system, figure [TT] 

In practice, we would not expect such precise stability. 
The attractor does present a favoured configuration in one 
parameter space but, as we shall demonstrate late, Jeans' 
equation means that a variety of density profiles and dis- 
persion will fit it. Since the analysis performed to produce 
these plots uses bins of fixed radius the subtle changes in 
the density profile are a source of noise that is not present 
in the control plot. In practice, one uses these plots to iden- 
tify when the changes in anisotropy are roughly constant 
at a given radius (since the rate of convergence is not con- 
stant normally) and when the fluctuations in anisotropy are 
reasonably close to those in the control plot. 

Noting that the colour scales of the two plots are differ- 
ent, we see that the system is actually very stably converged. 
Fluctuations in the outer edges of the system are very close 
to the control plot while the central regions are a little more 
unstable inside of 0.01 kpc, with moments when the system 
is as stable as the control and moments when it is not. This 
is taken as reasonable indication that the system has evolved 
as far as it is going to. Similar plots for our other systems 
demonstrate that they are also fully converged. 

As noted previously there is a lot of freedom in design- 
ing perturbation algorithms with particular freedom allowed 
in the numerical ranges of certain critical parameters such 
as the range of the scaling factor and the length of time al- 
lowed for the flow phase. In order to examine the effect of 
parameter choices on the progression of our simulations we 
ran four simulations using different combinations of scaling 
factors and flow times and then stopped them after only a 
few perturbations. We can then see which simulations have 
evolved further in that time and which are slower. We stop 
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Figure 12. The changes in anisotropy over 30 kicks of the New- 
tonian, isotropic system. Note the difference in scale from figure 
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Figure 14. The anisotropy profile of a standard isotropic system 
during the initial conditions, after the kick has been applied but 
before the system has evolved and then subsequent profiles as the 
system relaxes. 
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Figure 13. Attractor plots for four different combinations of 
algorithm parameters: Using scaling factors of 0.5 < / < 1.5 with 
1 (black) and 3 (red) dynamical times for the flow and 0.0 < f < 
2.0 with 1 (blue) and 3 (green) dynamical times. Each simulation 
has only been perturbed 4 times by this point. This number is 
arbitrary and was chosen as it presented the clearest plot. 

the simulations before they reach the attractor because at 
that point all the simulations should be lying on top of each 
other anyway and any differences in how they progressed to 
that point would be lost. 

Looking at figure [13] we can see why it is important to 
thoroughly examine the choice of algorithm as some param- 
eters have a significant impact while others do not, even if 
one might expect otherwise. For example, we show here that 
allowing the system to equilibriate for longer periods of time 
actually has almost no effect on the speed of convergence. 
This is slightly counter to expectation but amply demon- 
strates both the stability of all the intermediate solutions 
and the fact that the attractor does require the system to 
be driven to it to some degree but that that the final point 



of convergence is special compared to all the intermediate 
steps. In our scheme we see that the most important influ- 
encing factor is the size of the range of scaling factors applied 
with larger kicks promoting faster convergence. 



6 MECHANISM FOR CONVERGENCE 
6.1 One shock cycle in detail 

In order to understand the development of the convergent 
behaviour we investigate the changing state of the system 
during a single flow phase. We see several notable changes 
that progress over the course of the flow. The most interest- 
ing effect is the propagation of the anisotropy which seems to 
start in the inner regions and then subsequently travel out- 
wards as a wave of radial anisotropy which can be seen in 
figure [T4] This is because radial anisotropy denotes a popu- 
lation of particles with unusually high radial velocities which 
means that those particles will now be traveling towards the 
outer edges of the system. 

Note that the wave is still very clear even after 1 Tdyn- 
However, since figure [14] is plotted as a function of radius, 
there is actually very little mass between the wave and the 
edge of the system. This means that, if the system is left 
to evolve further, the profile as a function of mass will not 
change significantly. This is supported by figure ^3] where 
allowing the system to flow for a further 2 Tdyn has no dis- 
cernable impact on the system's evolution towards the at- 
tractor. Hence we can say that these systems are as close 
to equilibrium as makes no difference and can, for our pur- 
poses, be treated as fully evolved and suitable for another 
perturbation. 

The perturbation scheme itself does not instantaneously 
perturb the anisotropy profile to any significant degree but 
we see that within a very small amount of time we have a 
pronounced peak over a range of radii which then proceeds 
to migrate to the outer regions of the system. The peak 
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Figure 15. The anisotropy profile of the very radial IC during 
the initial conditions, after the kick has been applied but before 
the system has evolved and then subsequent profiles as the system 
relaxes. 



originates from the highest density area of the system which 
follows from the statistical nature of the cause of the effect. 

What this shows us is that the perturbation sets up a 
population of strongly radial particles but in such a way that 
they are not simply aligned specifically to be now travel- 
ing radially. If the latter were the case then the anisotropic 
signal would show up in the profile the instant after the 
perturbation algorithm completes whereas we see that some 
time, here at least 0.01 Tdyn, must pass before the anisotropy 
presents itself. This is supported by the aforementioned sym- 
metry of the perturbation algorithm in its treatment of ve- 
locity vectors and the fact that the algorithm does not know 
what the radial velocity vector of a given particle even as it 
is, as we showed in the previous section, coordinate- system 
independent. 

The cause for this radial preference is due to the statis- 
tics involved in putting a particle on a radial orbit. Consider 
a particle on a circular orbit around a centre of mass. If 
we perturb that particle using our scheme then, unless that 
particle is lucky enough to have it's velocity components all 
scaled by unity, it will either gain or lose energy and end up 
on a more radial orbit either now falling towards the centre 
of it's orbit or being ejected outwards. If we want to return 
our test particle onto a more circularised orbit we must now 
be extremely fortunate and apply just the right scaling to 
each velocity component so that not only does it end up with 
the correct total energy but it is now moving tangentially to 
the centre of it's new orbit. This is far less likely than it just 
being put on another elliptical orbit and thus contributing 
to the radial anisotropy. 

We do see different behaviour from systems that are ini- 
tially more radially anisotropic than the attractor. In figure 
[15] we still see that the perturbation itself does not change 
anything significantly but we are then missing the 'wave' of 
anisotropy visible in figure |14| This is because the mecha- 
nism here has more in common with standard radial orbit 
instability (ROI) and thus the most significant effect seen 
in the profile is the inward collapse of the very radial outer 



Figure 16. Virial equilibrium at the end of every flow phase for 
the isotropic Newtonian simulations. 



edges. The outward 'wave' of high energy particles produced 
by the perturbation is still there but the lack of an isotropic 
background prevents it from standing out. It's presence can 
be inferred from the increased size of the system. Since stud- 
ies of ROI show that radial populations are highly suscep- 
tible to this kind of collapse ( | Barnes et a l. 2009), what we 
observe in figure [15] is not unexpected. 

The radially anisotropic outer parts of the system 
quickly collapse radially inwards (the yellow line from figure 
15 ) as the kick has robbed them of energy. As they move in- 
wards they dramatically increase the radial anisotropy of 
the region they pass through until they settle into more 
circular orbits at smaller radii (blue and green lines). The 
more isotropic regions in our system behave the same as the 
system from figure [14] and thus we have a general increase 
radial anisotropy anyway. Overall, we find the nucleus be- 
comes more radial at the expense of the envelope thanks to 
the migration of the outer regions. We still see a very ra- 
dial population at the extreme edge of the system as the 
highly energetic, almost ejected 'wave' population, seen so 
clearly in figure [14] is still present, just overshadowed by the 
collapse of the existing radial population. 



6.2 Overall effect of the perturbation 

While the presence of an attractor in the parameter space 
seems quite clear it is not obvious what mechanism is at 
work to evolve the system towards that attractor. We can 
identify some effects, such as the radial population discussed 
in the previous paragraph, but it does not appear that we are 
driving the system significantly out of equilibrium with our 
perturbation scheme. To within reasonable limits the system 
never leaves virial equilibrium. We know, from the nature 
of our conservation laws, that the kinetic energy is rigidly 
conserved and the potential energy is never affected so the 
start of the flow is always virialised if the preceding flow 
phase ended virialised. Plotting the virial ratio at the end 
of every flow phase in figure [16] for the isotropic Newtonian 
case shows that the system is always in equilibrium. 

In addition to this, there is no obvious analytical path 
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to producing this kind of anisotropy. If our perturbation 
scheme treats each independent velocity component equally 
then an initially isotropic distribution should be transformed 
to another equivalent isotropic distribution. If we begin from 
the definition of r from equation [8] but re- written for clarity 
here as: 



v x x + v y y + v z z 



(9) 



What we are explicitly conserving is the energy in the 
system i.e. the square of the velocities: 



'?r= ( 



v x x + v y y + v z z\ 2 - x 2 - y 2 - z 2 r . xy 



+ ■ 



(10) 

Since we are in an orthogonal system the cross terms 
are all which leaves us with: 



o i o V , 9 



(ii) 



So, in a large enough sample size, the random scale fac- 
tors that we apply will sum to unity and we will not alter 
v 2 r - The only distribution that we do alter is the distribution 
of the particles in energy momentum space where the per- 
turbation acts to spread out the distribution with successive 
shocks. This can be imagined as, in each bin, the particles 
with the highest and lowest velocities might get scaled to 
even high and even lower velocities respectively which has 
the effect of broadening the distribution in that bin. 

Figure ^] shows how the kinetic energy distribution 
spreads out while the distribution in J becomes slightly more 
compacted towards 0. The latter effect is because of the in- 
creased radial anisotropy; for a given energy in the final 
state, more of that energy is in the radial modes leading 
to a lower characteristic angular momentum for that en- 
ergy. Note the clear boundary marking an excluded zone in 
the upper-right quadrant. This arises due to two limiting 
effects, both aspects of the restriction on particles leaving 
the system. Energy is restricted simply in relation to the 
local binding potential which in turn is related to radius 
which can be plausibly proxied by angular momentum. An- 
gular momentum is restricted because, as touched on before, 
angular momentum for a given energy is limited by the mo- 
mentum of a circularised particle with that energy. These 
two effects create a decaying cut-off defining an allowed re- 
gion of values. 

Another behaviour that seems characteristic of the per- 
turbation is the splitting of the system into two populations 
of particles: a 'nucleus' of particles who move towards the 
inner regions of the system and a small population of 'en- 
velope' particles which move further out. This is most obvi- 
ous in the larger, more extended systems as the collapse of 
the larger, 'nucleus' population is very noticeable whereas in 
most simulations it is the 'envelope' population that makes 
its presence most apparent as the increase in physical size 
increases the dynamical timescale and thus the time the sim- 
ulation must run for. 
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Figure 17. Evolution of Eki n ~J number density contours for the 
isotropic Newtonian. The hexagonal effect in low density areas is 
an artifact of the binning process used for the contours. 



6.3 Path to convergence 

Since we now have a robust grasp of the phenomenology of 
the attractor the final step before attempting to understand 
the mechanisms at work is to summarise the convergence 
towards to the attractor by examining the flow of individual 
bins. First we see if the similarities in large scale behaviour 
between different simulations are mirrored in similarities at 
the scale of a single bin. For this section we use only those 
simulations carried out in Newtonian gravity. Given that 
the binning used by our analysis codes is radius-based not 
mass-based the necessary differences in density profile due to 
differing gravitational paradigms would cause a systematic 
error in any side-by-sid e co mparison. 

As we see in Figure 19 the 20 th radial bin from the cho- 
sen simulations all follow a very similar path over the course 
of their evolution towards the attractor. That we have cho- 



sen the 20* bin for this is arbitrary and largely for aesthetic 
reasons as there are 80 similar plots for the other bins all 
of which show similar agreements. The poorest agreement is 
found in the furthest bins as the low number of particles at 
those radii makes them more susceptible to statistics while 
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Figure 18. The development of the nucleus/envelope arrange- 
ment caused by the perturbation after 30 cycles. Note that the 
kink at The reason that the nucleus/envelope arrangement is pro- 
duced is related to the aforementioned preferential selection of 
radial orbits. The cause will be examined in detail in section f6. 41 




Figure 19. The position of the 20 th bin from the 6 Newtonian 
simulations thus far. It is not important which line represents 
which simulation as all that the plot needs to convey is that the 
bins all move along the same path. 



the agreement of the innermost bins is harder to demon- 
strate visually as they move only a small distance. 

With this information we can now construct a plot that 
shows how any radial bin should move if dropped into a 
simulation anywhere in our parameter space. We take all 
the bins in all simulations and put them in the same space 
before using a 2D histogram to bin that space. We create 
a simple 'velocity' of each bin according to it's change in 
position between one step and the next before placing those 
vectors in our histogram bins and averaging them together. 
We can now drop test particles in our 'bin velocity field' and 
see how they move. 

Figure [20] can be interpreted as showing the vector field 



Figure 20. A representation of the attractor in terms of a vector 
field. Hotter colours represent larger rates of change of beta, ^j, 
in the parameter space between perturbations. 



that the attractor creates in the parameter space and thus 
shows the paths that bins will move down if they are placed 
in a simulation. The plot nicely summarises both the ap- 
parent universality of the attractor and also allows for pre- 
diction of the state of intermediate, partially converged sys- 
tems. Note that the need for discretised particles and a lack 
of data about certain tangential anisotropics (as the simula- 
tions evolved through the space too fast) has created some 
blank areas around ft = —0.5,7 + K — ~ 5. This is unavoid- 
able and is due to the compromise between a scarcity of very 
tangentially anisotropic data points and a very large number 
of radially anisotropic points. 



6.4 Exploration of a possible mechanism - 
bimodal perturbation 

So far we have several clues as to what causes the convergent 
behaviour. We know that the characteristics of the initial 
system and the environment are irrelevant and that alter- 
ing the constants used in the perturbation only changes the 
intermediate steps and not the final conclusion. Using the 
fact that, as seen in section [5^2] changing the amount of time 
for the flow phase does not have any impact on the rate of 
convergence we can suppose that maybe the flow phase only 
acts to allow the system to find the nearest Jeans'-stable 
state and does not actually contribute to the attractor it- 
self beyond the small changes seen in section |6.1| We might 
then ask if the attractor is caused primarily by the repeated 
action of the perturbation scheme and that the attractor 
represents the closest Jeans'-stable state to some limiting 
point for the perturbation scheme. 

One potential insight into this comes from the realisa- 
tion that all of our algorithms are actually special cases of 
a broader class of perturbation that always gives rise to the 
same behaviour for well-understood reasons. First, we will 
run down the various permutations of the algorithm and 
show why they all accomplish the same thing and can be 
thought of as being related. After establishing that one par- 
ticular algorithm can be used to represent all of the other 
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variants tested so far we will show why we would expect that 
algorithm to lead to some convergent state. 

As we have previously discussed, the decision to imple- 
ment a perturbation that applies different factors on each 
velocity component has the principle effect of slowing down 
the convergence as the shock is not as strong as that pro- 
vided by a similar scheme that uses the same factor for all 
three. Apart from this difference in convergence rate, how- 
ever, the algorithms have the same effect and underlying 
principle and thus we can think of the former being a spe- 
cial case of the latter. 

With that said, imagine a fairly weak perturbation (of 
the latter variety) where 0.75 < / < 1.25 and consider the 
effect of applying this perturbation n times in succession 
whilst allowing no opportunity for relaxation in between. 
We can see that, statistically, a small population will feel a 
perturbation like 0.75 n < f < 1.25 n which will correspond 
to a more pronounced perturbation. We then posit a limit- 
ing case for these perturbations where n cycles of any given 
perturbation will have the impact of a single perturbation 
where < / < 2, the strongest perturbation that we have 
employed so far and strongest possible symmetric perturba- 
tion. 

So, all perturbations of the form 1 — C < f < 1 + C 
can be thought of as reducing to < / < 2 over a sufficient 
number of iterations. Now we apply the same line of think- 
ing; what do multiple applications of this algorithm reduce 
to? In this instance, a small population of particles will end 
up with velocities of almost zero while a small population 
will experience a significant boost. This time, when we ap- 
ply the algorithm multiple times we will not be significantly 
change the energy in the population at v ~ as the most we 
can do is double the particles' velocity but at worst we can 
set it to exactly zero. This means that every time we per- 
form the perturbation we increase the size of the low speed 
population while, in order to preserve energy conservation, 
we pump more and more energy into a small population of 
escapers. 

The conclusion of this line of reasoning is that, in an 
extreme, theoretical limit, all of the perturbation schemes 
that perturb by this kind of random scaling cause a large 
population go undergo radial-infall as their energy is slowly 
sapped away and fed into a smaller population of particles on 
high energy orbits that form an envelope. This interpretation 
is supported by previous evidence showing an increase in 
central density even as the outer edges of the bound system 
move further and further out. To test this in the most violent 
limit we implement an algorithm that implements a bimodal 
scaling rather than a uniform one; particles are either put 
into radial infall by having their speed set to almost zero or 
are boosted into the high energy envelope. An example of 
such a perturbation could be that 90% of the particles are 
induced to radial-infall but the 10% envelope has 10 times 
it's original energy in order to make up conservation. 

For our tests an algorithm was designed that either in- 
duced radial infall or set the particle traveling at it's escape 
speed in such a way that, in the statistical limit, energy 
conservation is implied. So that the radial infall is not too 
perfect and to allow for the effects of the pre-existing veloc- 
ity dispersion to have an effect we set the lower limit to be 
just more than zero; s = 0.01. 




Figure 21. The almost immediate convergence using the bi- 
modal algorithm showing the initial isotropic system (black), the 
system after only one perturbation (blue) and after two pertur- 
bations (green) 



where / = 



Void 



' Void X 
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if rand u [0, 1] > / 
if rand u [0, 1] < / 



(12) 



(13) 



What we see in figure [2l] is striking. After only one 
perturbation the system is almost exactly on top of the at- 
tractor and only one perturbation later the system has com- 
pletely converged. It is interesting that the system appears 
to 'overshoot' the attractor after the first perturbation and 
only settle onto it after the second kick. This is possibly 
because the initial perturbation is so unrealistically violent 
that it simply takes longer than the allowed 3 Td yn to equi- 
libriate. 

This suggests that the prior line of reason is valid in so 
far as all perturbation algorithms that we have dealt with so 
far are actually special, slower cases of the bimodal scheme 
and that the attractor can be generated simply by forcing 
the system to undergo our radial-infall-plus-ejection scheme 
and then letting the system become Jeans' stable again. 

This gives us a natural explanation for why the initial 
anisotropy profile and theory of gravity do not affect the at- 
tractor under such perturbation schemes. Initial anisotropy 
does not matter as the perturbation loses that information 
with every successive perturbation by, in the above theoret- 
ical limit, putting all the particles on highly radial orbits 
by either catapulting them out to the fringes of the system 
or by instantaneously stopping all motion. Gravity makes 
no difference as the attractor here is formed by radial infall 
which is not a behaviour who's existence is dependent on 
the specifics of the theory of gravity beyond the details of 
the speeds involved. 

This also explains interesting results from HJS where 
a perturbation that acts only on the radial velocity com- 
ponents will not lead to convergence whereas one that acts 
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Before continuing it is important to point out that the 
mechanism we use is not pure ROI as it forces the system 
to undergo an ROI-like collapse regardless of the initial sus- 
ceptibility of the system to ROI. For example, simulations 
in Bellovary et al. (2008) that used isotropic initial con- 



ditions found them to be stable against ROI compared to 
radially anisotropic ones and Barnes et al. (2009) found a 



strong tendency towards the development of triaxiality. Our 
simulations evolve regardless of initial anisotropy and sec- 
tion [4?2] showed our simulations remain generally spherical 
as our perturbation has no preferred axis. The triaxiality in 
the deep MOND simulations is interesting but ultimately 
a statistical effect of having a much more extended system 
to start with. Accordingly, the following are currently pre- 
sented as interesting connections between our work and ROI 
with the strength of the connection a topic of ongoing in- 
vestigation. 



Figure 22. Altered shape of the attractor due to a strong radial 
bias in the perturbation algorithm. 



only on the tangential components will. In the tangential- 
only case, we can take the limiting cases where the particle 
ends up with no angular momentum and ends up on a purely 
radial orbit or is boosted onto a highly elliptical orbit. This 
picture is very similar to the infall mechanism and thus the 
same outcome is to be expected. However, if only the radial 
components are changed then the angular momentum of the 
particle cannot be removed. This means we have the same 
high energy limiting case but now the low energy limiting 
case is a circular orbit. This leads to a very radial outer en- 
velope, as before, but a distinct tangential anisotropy in the 
nucleus regions that can be seen as a distinct 'S-shape' bend 
in the anisotropy curve. 

This kind of behaviour can emerge in any code with a 
strong radial bias in the perturbation. The results in figure 
|22| come from runs which attempted to conserve both an- 
gular momentum and kinetic energy. After conserving an- 
gular momentum the only way to conserve energy was to 
only alter the radial velocity components of each particle so 
as not to disturb the momentum as simple iteration towards 
some unrestricted and fully self-consistent solution could not 
guarantee convergence and, if it did, would have lost a large 
amount of information about the form of the intended per- 
turbation. 



6.5 Connection to Radial Orbit Instability (ROI) 

In the previous section we have shown that the perturba- 
tion schemes being used in our work all appear to reduce 
the problem to one of the evolution of a system under ROI. 
Much work has already been done in studying whether ROI 
is the explanation for the universality of the NFW profile in 
dark matter halos (Macmillan & Henriksen 2006; Bellovary 



et al. 2008; Lapi & Cavaliere 2011) and several universal pro- 
files and convergent behaviours have emerged. While these 
results are not exclusive to ROI they are all related to the 
mechanisms of collapse and, given the discussion of the pre- 
vious section, a very radial mode of collapse, such as ROI, 
would seem to be the most promising lens through which to 
view our results. 



6.5.1 Universal beta profile 

In many simulations it was observed that the anisotropy 
profile would tend towards a common profile as well, al- 
though more complex than a clean power-law (Bellovary 



et al. 120081 iLapi & Cavaliere|2011| , with an isotropic nucleus 



that rises smoothly to radial anisotropy in the outer edges. 
There is some tolerance within this, however, as Bellovary] 
et al. (2008) found that the initial anisotropy profile did 



leave a lingering impression on the final profile 

We so find good convergence among the Newtonian sim- 
ulations although again we find that the MOND simulations 
stand out somewhat in figure [23] as well as the extremely ra- 
dial simulations. We would expect this kind of convergence 
regardless of any discussion about ROI as we have already 
discussed in detail the physics that leads to these profiles. 



6.5.2 Convergence of p /a 3 

One important result, although not one exclusive to ROI, 
was that systems formed by mergers (such as DM halos) dis- 
play convergence in a quantity that proxies phase-space den- 
sity ( [Taylor fe Navarro||2001| |Dehnen fe McLaughiin]|2005l ): 



(14) 



As noted in the references, the peculiar thing about this 
convergence is that neither p nor a 3 are themselves conver- 
gent, to a power-law or otherwise. However, the significance 
of this quantity does suggest that there may be some con- 
nection between it and our 7 + k axis. Here, in figure [24] all 
we investigate is whether our simulations develop the same 
power-law during convergence. 

Figure [24] shows us firstly that all the models start of off 
with very similar phase-space density profiles with the ex- 
ception of the extremely radial case. However, all the New- 
tonian models, after evolution, tend towards the r~ 2 (equa- 
tionfl4| power-law reported by Taylor &; Navarro (2001) and 



Dehnen & McLaughlin ( 2005 ) and in fact end up with pro- 
files that are closer together than their already good agree- 
ment to begin with. Standing out from this close conver- 
gence is the extremely radial case which ends up with a 
translated r~ 2 power-law. Also interesting is that this tight 
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Figure 23. Evolution of the anisotropy profile in the Newtonian 
isotropic (black), radial (blue), tangential (green) and extremely 
radial (orange) models as well as the MONDian isotropic (purple) 
and radial (blue) models from their initial states (dotted lines/top 
panel) to their overall final states (solid lines/bottom panel). 



Figure 24. Evolution of phase-space density in the Newtonian 
isotropic (black), radial (blue), tangential (green) and extremely 
radial (orange) models as well as the MONDian isotropic (purple) 
and radial (blue) models from their initial states (dotted lines/top 
panel) to their overall final states (solid lines/bottom panel). 



convergence is completely absent from the MONDian pro- 
files which retain good agreement with each other but have 
barely evolved from their initial conditions. 

What this shows is that our perturbation, developed 
from HJS's attempt to simply model halo mergers, does 
seem to display some characteristics of a system that has 
undergone repeated mergers. Accordingly, we may use this 
power-law as a possible avenue of investigation into our at- 
tractor formula. Particularly in light of the notable connec- 
tions to the next result. 



6.5.3 /3 - 7 relation 

Seeing as though we now have two convergent properties, 
both involving velocity dispersions and density, it is not 
surprising that there is a relation between both f3 and 7, 
the density slope, that is also a convergent result. Many 
studies have found evidence for a relationship between the 
two |Huss et al. 1 1999 



Bellovary et al. 2008; Lapi & Cavaliere 2011) and [Hansen 
&; Moore (2006|) went as far as positing a relationship of 



|2006| |Hansen et al.| 



Barnes et al.[2005| |Hansen Moore 
2006| IMacmillan & Henriksen| |200^ 



f3(r) — —0.15 + O.27. Although this result was not linked 
directly to ROI, instead being found in every relaxed struc- 
ture, it should still apply and several papers have used it for 
such. 

So, as we can see there are already several 'attractor' 
solutions relating to 7 and f3 in various parameter spaces if 
a system undergoes ROI. With this in mind, it is perhaps 
not surprising that convergent results are also a feature of 
our perturbations. We cannot draw too many conclusions as 
we know that our perturbations are not producing textbook 
ROI and, interestingly, our MONDian simulations seem to 
produce the attractor despite not adhering to various conver- 
gent behaviours normally seen in ROI. However, the amount 
of established relationships in ROI systems does lead to- 
wards the possibility of finding an analytical solution to the 
attractor from first-principles, at least in Newtonian grav- 
ity, as well as offering a clear set of criteria to test for the 
emergence of the attractor. 
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6.6 Ongoing tests 

Since we find that all the tests so far are related to this ra- 
dial infall mechanism to some degree we are designing tests 
that specifically exclude radial infall scenarios by not us- 
ing randomised rescaling of velocity components. Work is 
ongoing into a perturbation scheme where a velocity distri- 
bution function is chosen a priori and the system is forced 
to follow that curve. It meets the requirements of a pertur- 
bation in this context as it moves particles in velocity space 
into a non-equilibrium solution however, unlike previous al- 
gorithms, information about anisotropy should be preserved 
as only the total speeds are being effected and the predeter- 
mined distribution function should prevent the formation 
of the 'nucleus' and 'envelope' populations that drive the 
attract or in all previous examples. 

Care must be taken to prescribe a distribution function 
with well understood behaviour otherwise the behaviour of 
the attractor could be confused with eccentric behaviour 
caused by choosing an unphysical distribution function. 

Investigations are ongoing as to whether the results 
from ROI work can be used to predict an analytic solution 
to the attractor by narrowing down the relationship between 
our perturbation simulations, MOND and the more general 
body of work surrounding ROI phenomenology. While our 
systems are not exactly inducing ROI, there appear to be a 
large enough number of similarities for ROI to be a worth- 
while avenue of inquiry. 

A different set of perturbations are considered in Sparre 
& Hansen (2012, [in prep]), where the importance of violent 
relaxation is discussed. 



7 CONCLUSIONS 

We have tested a variety of initial conditions, gravity mod- 
els and perturbation schemes and found that in each case 
the system will tend towards a single set of solutions as de- 
scribed by HJS. The attractor appears independent of ini- 
tial anisotropy and prevailing theory of gravity due to the 
perturbations tested thus far all creating a system heavily 
influenced by radial-infall behaviours. This appears to be 
the only identifiable driving mechanism at this time as viri- 
alisation, Jeans' stability and the Antonov stability laws do 
not suggest any pressure to find new equilibria. 

As noted by HJS it is possible to create algorithms that 
defy the attractor such as the radial-only algorithms that 
have been mentioned previously. That particular case is not 
a cause for concern as it is clearly a very unrealistic per- 
turbation method. Work is ongoing into finding a realistic 
algorithm that, according to the infall approach, should not 
lead to convergence and to that end an algorithm based on 
prescribed velocity distribution functions is being designed. 

Finally, it is interesting to look at our density profiles 
compared against NFW profiles. This is purely a curiosity as, 
while it would be very interesting to find that the attractor 
is also an attractor in density space, there is no reason to 
expect it to be. First, it is important to note that, due to 
the different shape of the potential-density pair for MOND, 
the MONDian Plummer spheres have the same shape but a 
lower central density so that more matter can be present in 
the outer regions. 
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Figure 25. Final density profiles from a selection of simula- 
tions contrasted with an arbitrary NFW profile, arbitrary Hern- 
quist profile and the initial conditions of the Newtonian and weak 
MOND Plummer profiles. 



From figure [25] we see that the models have changed 
noticeably from their initial configurations. We see that the 
attractor does manifest some effect as the 'nucleus/envelope' 
arrangement is visible here as the increased density in the 
very centre and the outer regions has the effect of producing 
an almost Hernquist-like profile. The significantly steeper 
and higher central profiles for the Newtonian models com- 
pared to MOND are due to the higher core density of such 
models as, for instance, the very extended deep MOND pro- 
file has an apparently very small core density for much the 
same reason. The difference is due to the perturbation in- 
creasing the core density relative to it's initial state which 
is proportionally smaller in models with more pronounced 
MONDian effects. It is unlikely that this represents a state- 
ment about the behaviour of particles in a MONDian regime. 

Overall, this is an interesting and robust phenomena 
that we will continue to study. In the event that subsequent 
work demonstrates that the radial infall mechanism is only 
a small effect relevant to this class of perturbation then the 
consequences and opportunities for observational prediction 
are significant. 
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straightforward Newtonian case. There are a variety of dif- 
ferent functions that can be used, but all use x, the ratio 
of the strength of the acceleration field to the MONDian 
threshold value of ao = 3600(km/s) 2 /kpc, and have the 
same asymptotic behaviour: 

Put another way, for a given potential we can calcu- 
late the acceleration under Newtonian dynamics, gN, using 
Poisson's equation. We can relate this to the MONDian ac- 
celeration, g, using p(x). Here we use a common, simple 
p(x) (noting that we are still only concerned with spheri- 
cally symmetric problems): 

HL = Jf) =li(x)= ' (A3) 
g \a J 1 + x 

In order to generate theoretical acceleration curves, 
such as those seen in figure [6] a similar process is carried 
out. For the mu function described above, the process would 
be: 

^ = ^ = v{y) = 4-^ = - + S{\ + -) (A4) 

g N fJL(x) Ky) \a J 2^ VV 4 V ; 

Thus, any Newtonian acceleration profile can be related 
to a MONDian one by g = gNv(y)- In the case of the deep 
MOND systems a different p(x) is used as the system is 
entirely in the MOND regime. The deep MOND theoretical 
acceleration profiles seen in figure [6] are found as follows: 

V • [|V0|V0] = 47rGa p(r) (A5) 

(V0) 2 = 47rGa \ f p(r)r 2 dr (A6) 
r Jo 

V<t>=]/ GM{ X )a ° (A7) 

where, in our case, p(r) would be the profile for a Plum- 
mer model. 



APPENDIX A: 
DYNAMICS 



DETAILS OF MONDIAN 



For a spherical system the equation that NMODY uses ( Lon- 
drill o &; Nipoti|2008 ) to solve for gravitational attraction is 
flBekenstein fe Milgrom||1984| ): 



fi ( ^ J V0(r) 



AnGp(r) 



(Al) 



to give an acceleration field, g, for a potential <j>{r) re- 
sulting from a density profile p(r). The p(x) function is what 
dictates the strength of the MONDian modification to the 



APPENDIX B: METHODS FOR L 
CONSERVATION 

If we want to conserve angular momentum then, for exam- 
ple, we can require the Lx, Ly and Lz components of each 
bin to be equal before and after. The change in each mo- 
mentum component from before and after the perturbation 
is assessed and each particle given an equal share of that 
change, SL[x, y, z], to make up. Dealing in Cartesian coordi- 
nates, each velocity component affects two momentum com- 
ponents which means that each velocity component has two 
pieces of information about how it must change to move 
towards conservation. The subsequent change made to the 
velocity component is an average of these two values. For 
example, for v x : 
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5Ly\z\ 



+ 



-SLz\y\ 



(\x\ + \z\)z (\x\ + \y\)y 



(Bl) 



with cyclical permutations for the other components. 
It can be shown that this approach rapidly converges to the 
correct global value. A qualitative explanation is that, in ev- 
ery circumstance, the two velocity components furthest from 
their required values will move to improve the two worst an- 
gular momentum offsets. 



APPENDIX C: IMPACT OF SPHERICITY ON 
ANISOTROPY 

A simple thought experiment involves a flat disc with a pop- 
ulation of particles at all radii in circular orbits. If we then 
take that system and stretch it along one axis by perfectly 
setting up an injection of radial velocity for every particle 
at the correct time we will end up with a set of elliptical 
orbits and a non-circular system. If we now place a circular 
mask over the system for a binning procedure and look at 
the orbits of the particles we will find a mix of radial and 
tangential motion. The problem is that this will now look 
like a circular system with a degree of anisotropy rather than 
like a very regular elliptical system. 

This paper has been typeset from a TgX/ ETgX file prepared 
by the author. 



